% Lab 3 data processing

close all;
clear all;
clc;

angles = [-15 15 25 45 75 110];

info = csvread('relevant_sources.csv');
d65 = info(1:10:end, 2);
xbar = info(1:10:end, 3);
ybar = info(1:10:end, 4);
zbar = info(1:10:end, 5);
xbarsmooth = info(:, 3);
ybarsmooth = info(:, 4);
zbarsmooth = info(:, 5);
clear info;

data = csvread('lab_3_data_pruned.csv')';
[~, num_samps] = size(data);
wavelengths = data(:,1);
data(:,1) = ones(1, numel(wavelengths))'*100;
data = data/100; % Scale to proportion

% Calculate stuff
normalize = 100/sum(d65.*ybar);
X_vals = sum(data .* repmat(d65, [1, num_samps]) .* repmat(xbar, [1, num_samps]), 1).*normalize;
Y_vals = sum(data .* repmat(d65, [1, num_samps]) .* repmat(ybar, [1, num_samps]), 1).*normalize;
Z_vals = sum(data .* repmat(d65, [1, num_samps]) .* repmat(zbar, [1, num_samps]), 1).*normalize;
x_vals = X_vals./(X_vals+Y_vals+Z_vals);
y_vals = Y_vals./(X_vals+Y_vals+Z_vals);

% Grab the source coordinates
x_source = x_vals(1);
y_source = y_vals(1);

X_source = X_vals(1);
Y_source = Y_vals(1);
Z_source = Z_vals(1);

% Prune off the light source from the data
data = data(:,2:end);
X_vals = X_vals(2:end);
Y_vals = Y_vals(2:end);
Z_vals = Z_vals(2:end);
x_vals = x_vals(2:end);
y_vals = y_vals(2:end);

%Calculating sRGB:
sRGBMat = [3.2404542 -1.5371385 -0.4985314;...
        -0.9692660  1.8760108  0.0415560;...
        0.0556434 -0.2040259  1.0572252];
sRGB_vals = (sRGBMat * [X_vals/100;Y_vals/100;Z_vals/100])/1.4;


% Calculating LAB
X_norm = X_vals./X_source;
Y_norm = Y_vals./Y_source;
Z_norm = Z_vals./Z_source;

L = zeros(size(X_norm));
L(Y_norm <= 0.008856) = 903.3*Y_norm(Y_norm <= 0.008856);
L(Y_norm > 0.008856) = 116*Y_norm(Y_norm > 0.008856).^(1/3)-16;

% VALUE OF a*
a = 500 * (f(X_norm) - f(Y_norm));

% VALUE OF b*
b = 200 * (f(Y_norm) - f(Z_norm));




% Calculate the spectrum locus
spectrum_locus_x = xbarsmooth./(xbarsmooth+ybarsmooth+zbarsmooth);
spectrum_locus_y = ybarsmooth./(xbarsmooth+ybarsmooth+zbarsmooth);
spectrum_locus_x(end+1) = spectrum_locus_x(1);
spectrum_locus_y(end+1) = spectrum_locus_y(1);

% Draw Camson Chromaticities
figure;
hold on;
grid on;
axis square;
box on;
set(gca, 'XTick', 0:0.1:1);
set(gca, 'YTick', 0:0.1:1);
plot(spectrum_locus_x, spectrum_locus_y, 'k', 'LineWidth', 2);
plot(x_source, y_source, 'ok');
% Plot the blue chromaticities    
plot(x_vals(1:6), y_vals(1:6), '^b');
plot(x_vals(7:12), y_vals(7:12), '^g');
plot(x_vals(13:18), y_vals(13:18), '^r');

plot(x_vals(19:24), y_vals(19:24), 'xb');
plot(x_vals(25:30), y_vals(25:30), 'xg');
plot(x_vals(31:36), y_vals(31:36), 'xr');

plot(x_vals(37:42), y_vals(37:42), 'ob');
plot(x_vals(43:48), y_vals(43:48), 'og');
plot(x_vals(49:54), y_vals(49:54), 'or');

xlabel('x Chromaticity');
ylabel('y Chromaticiy');
title('x-y Chromaticities for Camson Paper Sets');
print('-depsc', 'paper_xy.eps');
% Note trend increases from sample 1 to 6 for all paper sets.










% Draw Camson a*b*
figure;
hold on;
grid on;
box on;
set(gca, 'XTick', -60:10:60);
set(gca, 'YTick', -40:10:40);
axis equal;
axis([ -60 60 -40 40])
plot(a(1:6), b(1:6), '^b');
plot(a(7:12), b(7:12), '^g');
plot(a(13:18), b(13:18), '^r');

plot(a(19:24), b(19:24), 'xb');
plot(a(25:30), b(25:30), 'xg');
plot(a(31:36), b(31:36), 'xr');

plot(a(37:42), b(37:42), 'ob');
plot(a(43:48), b(43:48), 'og');
plot(a(49:54), b(49:54), 'or');

xlabel('a*');
ylabel('b*');
title('L*a*b* Chromaticities for Camson Paper Sets');
% Note trend increases from sample 1 to 6 for all paper sets.
print('-depsc', 'paper_lab.eps');


% 
figure;
hold on;
for i = 1:6:72
    plot(angles, L(i:i+5), '-x', 'Color', sRGB_vals(:,i));
end
xlabel('Measurement Angle');
ylabel('L*');
title('L* Trends by Sample');
print('-depsc', 'l_trend.eps');




% Black xy
figure;
hold on;
grid on;
axis square;
box on;
set(gca, 'XTick', 0:0.1:1);
set(gca, 'YTick', 0:0.1:1);
plot(spectrum_locus_x, spectrum_locus_y, 'k', 'LineWidth', 2);
plot(x_source, y_source, 'ok');
% Plot the blue chromaticities    
plot(x_vals(55:60), y_vals(55:60), '^k');
plot(x_vals(61:66), y_vals(61:66), 'xk');
plot(x_vals(67:72), y_vals(67:72), 'ok');
% No trend for black patches.
print('-depsc', 'metal_xy.eps');






figure;
hold on;
Colors = {'r', 'g', 'b', 'c', 'm', 'y'};
for i = (1:6) + 0 
    plot(wavelengths, data(:,i), 'Color', sRGB_vals(:,i));
    plot(wavelengths(1:10:end), data(1:10:end,i), '^', 'Color', sRGB_vals(:,i), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
    plot(wavelengths, data(:,i+18), 'Color', sRGB_vals(:,i+18));
    plot(wavelengths(3:10:end), data(3:10:end,i+18), 's', 'Color', sRGB_vals(:,i+18), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
    plot(wavelengths, data(:,i+36), 'Color', sRGB_vals(:,i+36));
    plot(wavelengths(6:10:end), data(6:10:end,i+36), 'o', 'Color', sRGB_vals(:,+36), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
end
xlabel('Wavelength (nm)');
ylabel('Reflectance Factor');
title('Blue Paper Spectra by Angle');
print('-depsc', 'blue_spectra.eps');






figure;
hold on;
Colors = {'r', 'g', 'b', 'c', 'm', 'y'};
for i = (1:6) + 6 
    plot(wavelengths, data(:,i), 'Color', sRGB_vals(:,i));
    plot(wavelengths(1:10:end), data(1:10:end,i), '^', 'Color', sRGB_vals(:,i), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
    plot(wavelengths, data(:,i+18), 'Color', sRGB_vals(:,i+18));
    plot(wavelengths(3:10:end), data(3:10:end,i+18), 's', 'Color', sRGB_vals(:,i+18), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
    plot(wavelengths, data(:,i+36), 'Color', sRGB_vals(:,i+36));
    plot(wavelengths(6:10:end), data(6:10:end,i+36), 'o', 'Color', sRGB_vals(:,+36), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
end
xlabel('Wavelength (nm)');
ylabel('Reflectance Factor');
title('Green Paper Spectra by Angle');
print('-depsc', 'green_spectra.eps');






figure;
hold on;
Colors = {'r', 'g', 'b', 'c', 'm', 'y'};
for i = (1:6) + 12 
    plot(wavelengths, data(:,i), 'Color', sRGB_vals(:,i));
    plot(wavelengths(1:10:end), data(1:10:end,i), '^', 'Color', sRGB_vals(:,i), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
    plot(wavelengths, data(:,i+18), 'Color', sRGB_vals(:,i+18));
    plot(wavelengths(3:10:end), data(3:10:end,i+18), 's', 'Color', sRGB_vals(:,i+18), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
    plot(wavelengths, data(:,i+36), 'Color', sRGB_vals(:,i+36));
    plot(wavelengths(6:10:end), data(6:10:end,i+36), 'o', 'Color', sRGB_vals(:,+36), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
end
xlabel('Wavelength (nm)');
ylabel('Reflectance Factor');
title('Red Paper Spectra by Angle');
print('-depsc', 'red_spectra.eps');






figure;
hold on;
Colors = {'r', 'g', 'b', 'c', 'm', 'y'};
for i = (1:6) + 54
    plot(wavelengths, data(:,i), 'Color', sRGB_vals(:,i));
    plot(wavelengths(1:10:end), data(1:10:end,i), '^', 'Color', sRGB_vals(:,i), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
    plot(wavelengths, data(:,i+6), 'Color', sRGB_vals(:,i+6));
    plot(wavelengths(3:10:end), data(3:10:end,i+6), 's', 'Color', sRGB_vals(:,i+6), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
    plot(wavelengths, data(:,i+12), 'Color', sRGB_vals(:,i+12));
    plot(wavelengths(6:10:end), data(6:10:end,i+12), 'o', 'Color', sRGB_vals(:,+12), 'MarkerFaceColor', Colors{mod(i-1,6)+1});
end
xlabel('Wavelength (nm)');
ylabel('Reflectance Factor');
title('Black Metal Spectra by Angle');
print('-depsc', 'metal_spectra.eps');
% Spectral plot for each angle for the three blacks

